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Abstract 

We report extensive Monte Carlo simulations of the Widom-Rowlinson lattice model 
in two and three dimensions. Our results yield precise values for the critical activities 
and densities, and clearly place the critical behavior in the Ising universality class. 
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The Widom-Rowlinson (WR) hard-sphere mixture is perhaps the simplest binary fluid 
model showing a continuous unmixing transition, and has been the subject of considerable 
study regarding its thermal and interfacial properties [p], 0, as has the Gaussian /-function 
version of the model introduced somewhat earlier by Helfand and Stillinger || . Despite this 
interest, however, definitive results on the location and nature of the WR critical point are 
lacking. Indeed, except for the essentially inconclusive series analyses of the Gaussian version 
0, ||, little has appeared by way of precise quantitative analysis on the critical behavior of 
any WR-type model. As a first step in this direction we have performed extensive simulations 
of the lattice-gas analog of the WR hard-sphere mixture the Widom-Rowlinson lattice model 
(WRL) @. 

In the original WR model, AB pairs interact via a hard-sphere potential whilst AA 
and BB pairs are noninteracting. By the WRL model we mean a two-component lattice 
gas in which sites may be at most singly occupied, and in which nearest-neighbor A-B 
pairs are forbidden. Like the WR model, this is evidently an athermal model (all allowed 
configurations are of the same energy), and is characterized solely by the densities of the two 
species or the corresponding chemical potentials \ia and \lb- (/m = Hb = A*c of course, at the 
critical point.) The WRL may be viewed as an extreme member of a family of binary alloy 
models. The Ising model, as is known, may be transcribed into such a model by identifying 
up and down spins with A and B particles, resp., yielding a "close-packed" alloy which 
unmixes at the Ising critical temperature. Allowing a small fraction of vacant sites results in 
a "dilute binary alloy" (DBA) with a somewhat depressed critical temperature; continuing 
the dilution process, one arrives at a model with T c = 0. This zero-temperature terminus 
of the DBA critical line is precisely the WRL critical point. One is then led to ask whether 
the entire line shares a common critical behavior, or whether its character changes at some 
point. Although the former expectation is clearly favored on the basis of universality, a 
careful examination of this question nevertheless appears worthwhile. One should also note 
that whilst in the WRL there is no temperature per se, ha+^b is a temperature-like variable, 
so that along the symmetry line /i^ = ^b = A*> the susceptibility scales as \ — (a* ~ A^) 7 , 
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the order parameter, pa — Pb — (p — Pc) 13 , (p > Pc), and so on ||. {h = pa — Pb plays the 
role of an external field.) 

Our simulations employ periodic square or cubic lattices of side L (L = 10, 20, 40, 80, 
and 160 in two dimensions; 8, 12, 16, 24, 40, and 64 in three), and are performed in the 
grand canonical ensemble, along the symmetry line, za = zb = z = e^' kT . The state of site 
i is conveniently represented by cr; = 1, 0, or -1, corresponding to i occupied by A, vacant, or 
occupied by B. Three kinds of moves are employed: 

(i) A "flip" in which, if o x ^ 0, <j\ — > —a\ or with equal likelihood, while if o"; = 0, &\ — > ±1 
with equal likelihood. Moves vacating a site are accepted with probability min[l, 1/z], those 
adding a new particle with probability min[l, z], provided no A-B nearest neighbor pairs are 
formed. Acceptance of changes in species is solely contingent on the latter condition. 

(ii) Exchange of the states of a pair of randomly chosen sites (anywhere in the system). 
Acceptance again depends on no A-B pair being generated. 

(iii) A "cluster flip" in which the entire particle cluster connected to site i is changed to the 
opposite species. Since clusters are bounded exclusively by vacant sites, such moves never 
generate A-B pairs and are always accepted. 

In processes (i) and (iii) site i is chosen at random; of course if o\ = there is no cluster 
flip. Cluster flips are of particular simplicity in WR models, and allow for rapid relaxation 
of any species imbalance. They are, however, time-consuming, particularly in large systems 
near or above the critical point, since an appreciable fraction of sites may belong to a single 
cluster. We found a reasonably efficient procedure was to attempt such moves about ten times 
per lattice update (lud). (A lud denotes L d moves of types (i) and (ii).) Block averages (over 
100 lud's) of mole fractions were then equal to within about in part in 10 3 . Our simulations 
furnish the following quantities as a function of z: the density p and its variance; the order 
parameter m = \p a — Pb\ and its variance; the probability distribution P(A) of the population 
excess A = Na — N B ; the reduced fourth cumulant, U(z, L) = 1— < m 4 > /3 < m 2 > [|J. 
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We used two independent methods to locate the critical activity z c : order-parameter 
histograms and fourth-cumulant crossings. The former method is based on the observation 
that -P(A) undergoes a clear change in character at a certain activity z c (L): for z < z c (L), 
-P(A) is peaked at zero, whilst above z c (L) it is bimodal; just at z c it has a broad plateau. 
Finite-size scaling ideas J7|, |8| imply that as L — > oo, 

z c (L) ^z. + AL- 1 ^. (1) 

In analyzing the data we adopt, as a working hypothesis, the Ising class //-value. Plotting 
the two-dimensional z c (L) data vs. L~ l yields a straight line with intercept z c = 2.0644, 
and a similar analysis in 3-d (using v = 0.629 0) gives z c = 0.7842. The widely-employed 
fourth-cumulant analysis uses the fact that U(z c , L) rapidly converges to a nontrivial value 
U* with increasing L [|. In our studies the crossings of the three largest L- values agree to 
within uncertainty. In 2d the crossing is at z c = 2.0636, U* = 0.607, while in 3d we find 
z c = 0.7840, U* = 0.478(5). Our final estimates for the critical activities are: 

z c = 2.0604(4) (2d), (2) 

z c = 0.7841(1) (3d). (3) 

The corresponding critical densities are p c = 0.618(1) (2d) and 0.3543(1) (3d). (The Bethe- 
Guggenheim approximation yields 0.429 and 0.2941, respectively, for the critical densities in 
2 and 3d 0) The hypothesis of Ising-like critical behavior finds significant support in the 
fact that the z c (L) plots are linear (this confirms the assumed value of v), and from the good 
agreement of our U* with the accepted values of 0.61 in 2d || and 0.47 in 3d ||. Further 
support appears when we examine the order parameter and susceptibility. 

In Figs. 1 and 2 we plot, for 2d and 3d, respectively, the scaled order parameter, rh = 
If! v m vs. the scaling variable t = L l ' v \z — z c \, using the Ising model values for the 
exponents. The data are seen to collapse nicely, and for 2d are fully consistent with f3 = 1/8, 
i.e., the slope of the scaling plot is ~ 1/8 for z > z c , and ~ —7/8 for z < z c . The slope of 
the 3d data is a bit low, a fit to the most linear portion yielding about 0.31, as compared 
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with the accepted value (3 = 0.326(4) 0. Finite size effects seem the most likely cause for 
this small discrepancy. The susceptibility is given by x = d < A > /dh = \var[A]. To 
study its scaling we plot % = L'" 1 ^ vs. t = L~ x l v \z — z m (L)\, where z m (L) is the activity 
at which x{ z i L) takes its maximum. In this case there is a reasonably good data collapse 
and the slopes of the scaling plots are in good agreement with the accepted values of 7, as 
seen in Figs. 3 and 4. In WR systems the analog of the specific heat is the compressibility, 
given by k = L d var[p]/ p 2 . We are unable to obtain k to good precision in the critical region 
(similar difficulty is encountered for the specific heat in Ising simulations), but can at least 
report that in two dimensions the maximum compressibility grows oc InL, as expected for 
models in the Ising class. Finally, it is of some interest to know the shape of the coexistence 
curves in two and three dimensions; these are shown in Fig. 5. 

In summary, we have performed extensive and detailed simulations of the Widom-Rowlinson 
lattice model in two and three dimensions, and find strong evidence of Ising-like critical be- 
havior. We believe that further studies of other aspects of the WRL are in order, for example, 
surface free energies and roughening, critical dynamics, and the influence of a driving field. 
Series analysis of critical behavior also appears feasible and worthwhile. Such investigations 
promise to shed new light on the Ising universality class, as well as on critical behavior in 
binary fluids. 
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ergy Sciences, Office of Energy Research, U.S. Department of Energy, and G.S. acknowledges 
the support of the National Science Foundation. 
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Figure Captions 



Fig. 1. Scaling plot of the order parameter, 2d. Open circles: L = 20; open diamonds: 
L = 40; filled circles: L = 80; squares: L = 160. The straight lines have slopes of 1/8 and 
-7/8. 

Fig. 2. Scaling plot of the order parameter, 3d. Open diamonds: L = 8; open squares: 
L = 12; open circles: L = 16; filled circles: L = 24; filled squares: L = 40. The line has 
slope 0.326. 

Fig. 3. Scaling plot of the susceptibility, 2d. Symbols as in Fig. 1. The slope of the line is 
-1.75. 

Fig. 4. Scaling plot of the susceptibility, 3d. Symbols as in Fig. 2. The slope of the line is 
-1.247. 

Fig. 5. Coexistence curves in 2d (broken line) and 3d (solid line). 
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